#include "ImageFea.h"

/*
int ExtractSiftFea(const char* str_src, const char* str_dst)
{
    if(NULL==str_src || NULL==str_dst)
        return -1;
        
    IplImage * ptImage = NULL;
	if (NULL == (ptImage = cvLoadImage(str_src, 0))) {
		//AfxMessageBox("Can't open image " + imgName + "!");
		printf("Can't open image %s!",str_src);
		return -2;
	}

	feature * fea;
	int size = 0;
	size = sift_features(ptImage, &fea);
	if (size == 0) {
		cvReleaseImage(&ptImage);
		return -3;
	}
	
	FILE*fp=fopen(str_dst,"w");
	if(NULL!=fp)
	{
        int i,j;
        int siftLen = fea[0].d;
        float ratio = 10.0;
        float curScale;
        fprintf(fp, "%d %d \n", size, siftLen + 4);
        for (i = 0; i < size; i++) {
            curScale = ratio * fea[i].scl;
		    fprintf(fp, "%f %f %f %f ", fea[i].x, fea[i].y, curScale, fea[i].ori);
		    for (j = 0; j < siftLen; j++) {
			    fprintf(fp, "%f ", fea[i].descr[j]);
		    }
		    fprintf(fp, "\n");
	    }
        fclose(fp);
    }
	
	//for(int i=0;i<size;i++){
    //    printf("%d,%d\n",size,i);
    //    if(NULL!=fea[i].feature_data)
    //        free(fea[i].feature_data);
    //}
	//delete[]fea;    
	free(fea);
	
	cvReleaseImage(&ptImage);    
    return size;
}
*/

int ExtractGrayHist(const char* imageFilename, const char* histoFilename
, const int nGridV, const int nGridH, const int nHistL)
{
    if(imageFilename==NULL || histoFilename==NULL || nGridV<1 || nGridH<1 || nHistL<1)
        return -1;
        
    
    IplImage*   ptImage=::cvLoadImage(imageFilename);
    if(NULL==ptImage)
        return -2;
    CvSize      imgSize=::cvGetSize(ptImage);
    IplImage*   ptGray=::cvCreateImage(imgSize,8,1);
    if(ptImage->nChannels==3)
        cvCvtColor(ptImage,ptGray,CV_RGB2GRAY);
    else
        cvCopy(ptImage,ptGray);
    
    float *fea=new float[nGridV*nGridH*nHistL];
    memset(fea,0,sizeof(float)*nGridV*nGridH*nHistL);
    
    int interV=imgSize.height/nGridV;
    int interH=imgSize.width/nGridH;
    int band=256/nHistL;
    // extract features
    unsigned char *ptrSrc;
    float *ptrFea;
    //int value;
    int iV,iH,iGridV,iGridH,iX;
	for(iV=0;iV<imgSize.height;iV++)
	{
        iGridV=iV/interV;
        if(iGridV>=nGridV)
            iGridV=nGridV-1;
		ptrSrc=(unsigned char*)cvPtr2D(ptGray,iV,0,0);
		iX=0;
		for(iGridH=0;iGridH<nGridH;iGridH++)
		{
            ptrFea=fea+((iGridV*nGridH)+iGridH)*nHistL;
            for(iH=0;iH<interH;iH++,iX++)
            {
                ptrFea[*ptrSrc++/band]++;
            }
            if(iGridH==nGridH-1 && iX<imgSize.width)//the last several
            {
                for(;iX<imgSize.width;iX++)
                    ptrFea[*ptrSrc++/band]++;
            }
        }
	}
	float sum;
	ptrFea=fea;
	for(iGridV=0;iGridV<nGridV;iGridV++)
	{
        for(iGridH=0;iGridH<nGridH;iGridH++)
        {
            sum=0;
            for(iX=0;iX<nHistL;iX++)
                sum+=ptrFea[iX];
            for(iX=0;iX<nHistL;iX++)
                ptrFea[iX]/=sum;
            ptrFea+=nHistL;
        }
    }
    
    ::cvReleaseImage(&ptImage);
    ::cvReleaseImage(&ptGray);
    
    FILE *afile=fopen(histoFilename,"wb");
    fwrite(fea,sizeof(float),nGridV*nGridH*nHistL,afile);
    fclose(afile);
    
    delete[]fea;
    
    return 0;
}
int ExtractGrayHist
(IplImage*   ptImage, float *fea
, const int nGridV, const int nGridH, const int nHistL)
{
    if( NULL==ptImage || NULL==fea || nGridV<1 || nGridH<1 || nHistL<1)
        return -1;        

    CvSize      imgSize=::cvGetSize(ptImage);
    IplImage*   ptGray=::cvCreateImage(imgSize,8,1);
    if(ptImage->nChannels==3)
        cvCvtColor(ptImage,ptGray,CV_RGB2GRAY);
    else
        cvCopy(ptImage,ptGray);
    
    memset(fea,0,sizeof(float)*nGridV*nGridH*nHistL);
    
    int interV=imgSize.height/nGridV;
    int interH=imgSize.width/nGridH;
    int band=256/nHistL;
    // extract features
    unsigned char *ptrSrc;
    float *ptrFea;
    int iV,iH,iGridV,iGridH,iX;
	for(iV=0;iV<imgSize.height;iV++)
	{
        iGridV=iV/interV;
        if(iGridV>=nGridV)
            iGridV=nGridV-1;
		ptrSrc=(unsigned char*)cvPtr2D(ptGray,iV,0,0);
		iX=0;
		for(iGridH=0;iGridH<nGridH;iGridH++)
		{
            ptrFea=fea+((iGridV*nGridH)+iGridH)*nHistL;
            for(iH=0;iH<interH;iH++,iX++)
            {
                ptrFea[*ptrSrc++/band]++;
            }
            if(iGridH==nGridH-1 && iX<imgSize.width)//the last several
            {
                for(;iX<imgSize.width;iX++)
                    ptrFea[*ptrSrc++/band]++;
            }
        }
	}
	float sum;
	ptrFea=fea;
	for(iGridV=0;iGridV<nGridV;iGridV++)
	{
        for(iGridH=0;iGridH<nGridH;iGridH++)
        {
            sum=0;
            for(iX=0;iX<nHistL;iX++)
                sum+=ptrFea[iX];
            for(iX=0;iX<nHistL;iX++)
                ptrFea[iX]/=sum;
            ptrFea+=nHistL;
        }
    }
    
    ::cvReleaseImage(&ptGray);    
    return 0;
}
int ExtractGrayHist(const char* imageFilename, float *fea
, const int nGridV, const int nGridH, const int nHistL)
{
    if( NULL==imageFilename || NULL==fea || nGridV<1 || nGridH<1 || nHistL<1)
        return -1;
        
    IplImage*   ptImage=::cvLoadImage(imageFilename);
    if(NULL==ptImage)
        return -2;
    ExtractGrayHist(ptImage,fea,nGridV,nGridH,nHistL);
    
    ::cvReleaseImage(&ptImage);
    return 0;
}
// mean color value for each gird
int ExtractColorFea(string imageFilename, string colorFilename, int nGridV, int nGridH)
{
    if(imageFilename.empty() || colorFilename.empty() || nGridV<1 || nGridH<1)
        return -1;
    
    IplImage*   ptImage=::cvLoadImage(imageFilename.data());
    if(NULL==ptImage)
        return -2;
    CvSize      imgSize=::cvGetSize(ptImage);
    IplImage*   ptHSV=::cvCreateImage(imgSize,8,3);
    IplImage*   ptH=::cvCreateImage(imgSize,8,1);
    cvCvtColor(ptImage,ptHSV,CV_RGB2HSV);
    cvSplit(ptHSV,ptH,NULL,NULL,NULL);
    
    float *fea=new float[nGridV * nGridH];
    memset(fea,0,sizeof(float)*nGridV * nGridH);
    
    int interV=imgSize.height/nGridV;
    int interH=imgSize.width/nGridH;
    // extract features
    unsigned char *ptrSrc;
    float *ptrFea;
    //int value;
    int iV,iH,iGridV,iGridH,iX;
	for(int iV=0;iV<imgSize.height;iV++)
	{
        iGridV=iV/interV;
        if(iGridV>=nGridV)
            iGridV=nGridV-1;
		ptrSrc=(unsigned char*)cvPtr2D(ptH,iV,0,0);
		iX=0;
		for(iGridH=0;iGridH<nGridH;iGridH++)
		{
            ptrFea=fea+(iGridV*nGridH)+iGridH;
            for(iH=0;iH<interH;iH++,iX++)
            {
                *ptrFea+=*ptrSrc++;
            }
            if(iGridH==nGridH-1 && iX<imgSize.width)//the last several
            {
                for(;iX<imgSize.width;iX++)
                    *ptrFea+=*ptrSrc++;
            }
        }
	}
	
	ptrFea=fea;
    for(iV=0;iV<nGridV-1;iV++)
    {
        for(iH=0;iH<nGridH-1;iH++)
            *ptrFea++/=(interV*interH);
        *ptrFea++/=(interV*(imgSize.width-interH*(nGridH-1)));
    }
    for(iH=0;iH<nGridH-1;iH++)
        *ptrFea++/=((imgSize.height-interV*(nGridV-1))*interH);
    *ptrFea++/=((imgSize.height-interV*(nGridV-1))*(imgSize.width-interH*(nGridH-1)));
	
	//for(int i=0;i<8*8;i++)printf("%f\n",fea[i]);getchar();
	
    FILE *afile=fopen(colorFilename.data(),"wb");
    fwrite(fea,sizeof(float),nGridV*nGridH,afile);
    fclose(afile);
	delete[] fea;
	
    ::cvReleaseImage(&ptImage);
    ::cvReleaseImage(&ptHSV);
    ::cvReleaseImage(&ptH);
    
    return 0;
}

Gabor * constructGaborsLT(int orientation, int frequency)
{
	int i,j;
	//string name="gk/";//gabor kernel
	//mkdir( name.data() );
	//char buffer1[64],buffer2[64];
	
	Gabor *gaborsLT=new Gabor[orientation*frequency];

	for(i=0;i<orientation;i++)
	{
		for(j=0;j<frequency;j++)
		{
			gaborsLT[i*frequency+j].Init(PI*i/orientation, frequency_select[j*2], Sigma, F);
			//real
			//cvSave((name+"r_"+itoa(i,buffer1,10)+"_"+itoa(j,buffer2,10)).data(),gaborsLT[i*frequency+j].get_matrix(REAL));
			//image
            //cvSave((name+"i_"+itoa(i,buffer1,10)+"_"+itoa(j,buffer2,10)).data(),gaborsLT[i*frequency+j].get_matrix(IMAG));
		}
	}

	return gaborsLT;
}
int constructGaborsLT(Gabor * gabors, const int frequency, const int orientation)
{
	int i,j;
	if(gabors==NULL)
	   return -1;

	for(i=0;i<frequency;i++)
	{
		for(j=0;j<orientation;j++)
		{
			gabors[i*orientation+j].Init(PI*i/orientation, frequency_select[j], Sigma, F);
		}
	}
	
	return 0;
}
Gabor* loadGaborsLT(int orientation, int frequency)
{
	int i,j;
	string name="gk/";//gabor kernel
	
	Gabor *gaborsLT=new Gabor[orientation*frequency];
	//CvMat *pMat=NULL;

	for(i=0;i<orientation;i++)
	{
		for(j=0;j<frequency;j++)
		{
		}
	}

	return gaborsLT;
}
int ExtractGaborFea
(string imageFilename, string gaborFilename
, int nGridV, int nGridH
, int nGaborO, int nGaborF)
{
    if(imageFilename.empty() || gaborFilename.empty() || nGridV<1 || nGridH<1 || nGaborO<1 || nGaborF<1)
        return -1;
    
    IplImage*   ptImage=::cvLoadImage(imageFilename.data());
    if(NULL==ptImage)
        return -2;
        
    CvSize      imgSize=::cvGetSize(ptImage);
    IplImage* currImage=NULL;
    if( imgSize.width*imgSize.height > 400*300 )
    {
        if(imgSize.width>imgSize.height)
            currImage=::cvCreateImage(cvSize(800,600),ptImage->depth,ptImage->nChannels);
        else
            currImage=::cvCreateImage(cvSize(600,800),ptImage->depth,ptImage->nChannels);        
        ::cvResize(ptImage,currImage,CV_INTER_AREA);
        
    }
    else  
        currImage=::cvCloneImage(ptImage);    
    imgSize=::cvGetSize(currImage);
    
    IplImage*   ptGray=::cvCreateImage(imgSize,8,1);
    cvCvtColor(currImage,ptGray,CV_RGB2GRAY);
    IplImage*   ptConv= cvCreateImage(imgSize,IPL_DEPTH_32F, 1);
    
    Gabor *gaborsLT;
    gaborsLT=constructGaborsLT(nGaborO, nGaborF);
    //gaborsLT=loadGaborsLT(nGaborO, nGaborF);

    float* fea=new float[nGaborO*nGaborF*nGridV*nGridH*2];
    float* ptrfea,*ptrfea0=NULL;
    memset(fea,0,sizeof(float)*nGaborO*nGaborF*nGridV*nGridH*2);
    int interV=imgSize.height/nGridV;
    int interH=imgSize.width/nGridH;
    int i,j,ix,iy,iV,iH;//,temp;
    unsigned char* ptr;
	for(i=0;i<nGaborO;i++)
	{
		for(j=0;j<nGaborF;j++)
		{
            gaborsLT[i*nGaborF+j].conv_img(ptGray, ptConv, REAL);
            
            ::cvNamedWindow("it");
            ::cvShowImage("it",ptConv);
            ::cvWaitKey();
            ::cvDestroyWindow("it");
            
            ptrfea0=fea+(i*nGaborF+j) *nGridV*nGridH*2;  
            //mean         
            for(iy=0;iy<imgSize.height;iy++)
            {
                iV=iy/interV;
                if(iV>=nGridV)
                    iV=nGridV-1;
                ptr=(unsigned char*)cvPtr2D(ptConv,iV,0,0);
                ptrfea=ptrfea0+iV*nGridH*2;
                for(ix=0;ix<imgSize.width;ix++)
                {
                    iH=ix/interH;
                    if(iH>=nGridH)
                        iH=nGridH-1;
                    ptrfea[iH+iH]+=*ptr++;
                }
            }
            for(int i=0;i<64;i++)printf("%f\n",fea[i]);getchar();
            ptrfea=ptrfea0;
            for(iV=0;iV<nGridV-1;iV++)
            {
                for(iH=0;iH<nGridH-1;iH++)
                {
                    *ptrfea/=(interV*interH);
                    ptrfea+=2;
                }
                *ptrfea/=(interV*(imgSize.width-interH*(nGridH-1)));
                ptrfea+=2;
            }
            for(iH=0;iH<nGridH-1;iH++)
            {
                *ptrfea/=((imgSize.height-interV*(nGridV-1))*interH);
                ptrfea+=2;
            }
            *ptrfea/=((imgSize.height-interV*(nGridV-1))*(imgSize.width-interH*(nGridH-1)));
            //variation
            for(iy=0;iy<imgSize.height;iy++)
            {
                iV=iy/interV;
                if(iV>=nGridV)
                    iV=nGridV-1;
                ptr=(unsigned char*)cvPtr2D(ptConv,iV,0,0);
                ptrfea=ptrfea0+iV*nGridH*2;
                for(ix=0;ix<imgSize.width;ix++)
                {
                    iH=ix/interH;
                    if(iH>=nGridH)
                        iH=nGridH-1;
                    ptrfea[iH+iH+1]+=((ptrfea[iH+iH]-*ptr)*(ptrfea[iH+iH]-*ptr));
                    ptr++;
                }
            }
            ptrfea=ptrfea0;
            for(iV=0;iV<nGridV-1;iV++)
            {
                for(iH=0;iH<nGridH-1;iH++)
                {
                    ptrfea[1]=sqrt(ptrfea[1]/(interV*interH));
                    ptrfea+=2;
                }
                ptrfea[1]=sqrt( ptrfea[1]/(interV*(imgSize.width-interH*(nGridH-1))) );
                ptrfea+=2;
            }
            for(iH=0;iH<nGridH-1;iH++)
            {
                ptrfea[1]/=((imgSize.height-interV*(nGridV-1))*interH);
                ptrfea+=2;
            }
            ptrfea[1]/=((imgSize.height-interV*(nGridV-1))*(imgSize.width-interH*(nGridH-1)));
            /**/
		}
	}
    FILE *afile=fopen(gaborFilename.data(),"wb");
    fwrite(fea,sizeof(float),nGaborO*nGaborF*nGridV*nGridH*2,afile);
    fclose(afile);
	delete[]fea;
	
    delete[] gaborsLT;
    ::cvReleaseImage(&ptImage);
    ::cvReleaseImage(&ptGray);
    ::cvReleaseImage(&ptConv);
    ::cvReleaseImage(&currImage);
    
    return 0;
}
int ExtractGaborFea(const char* imageFilename, const char* gaborFilename, Gabor* gabors
, const int nGaborF, const int nGaborO, const int nGridV, const int nGridH)
{
    if(imageFilename==NULL || gaborFilename==NULL || NULL==gabors || nGridV<1 || nGridH<1 || nGaborO<1 || nGaborF<1)
        return -1;
    
        
    IplImage*   ptImage=::cvLoadImage(imageFilename);
    if(NULL==ptImage)
        return -2;
    
    // resize to tiny image (64,64)
    IplImage* ptTinyImage=cvCreateImage(cvSize(64,64),8, ptImage->nChannels);
    cvResize( ptImage, ptTinyImage);
        
    // convert to gray image
    CvSize      imgSize=::cvGetSize(ptTinyImage);
    IplImage*   ptGray=::cvCreateImage(imgSize,8,1);
    if(ptTinyImage->nChannels==3)
        cvCvtColor(ptTinyImage,ptGray,CV_RGB2GRAY);
    else
        cvCopy(ptTinyImage,ptGray);
        
    IplImage*   ptConv= cvCreateImage(imgSize,IPL_DEPTH_32F, 1);
    
    float *fea=new float[nGaborF*nGaborO*nGridV*nGridH];
    memset(fea,0,sizeof(float)*nGaborF*nGaborO*nGridV*nGridH);
    
    float *ptrfea=NULL;
    float *ptrfea1=NULL;
    float* ptr=NULL;
    int interV=imgSize.height/nGridV;
    int interH=imgSize.width/nGridH;
    // for each frequency
    int iFre,iOri,iV,iH;
    for(iFre=0;iFre<nGaborF;iFre++)
    {
        // for each orientation
        for(iOri=0;iOri<nGaborO;iOri++)
        {
            ptrfea=fea+(iFre*nGaborO+iOri)*nGridV*nGridH;
            gabors[iFre*nGaborO+iOri].conv_img(ptGray, ptConv, MAG);
            
            /*
            ::cvNamedWindow("it");
            ::cvShowImage("it",ptConv);
            ::cvWaitKey();
            ::cvDestroyWindow("it");
            */
            
            //means         
            float sum=0;
            for(int iy=0;iy<imgSize.height;iy++)
            {
                iV=iy/interV;
                if(iV>=nGridV)
                    iV=nGridV-1;
                // for each row
                ptr=(float*)cvPtr2D(ptConv,iV,0,0);
                ptrfea1=ptrfea+iV*nGridH;
                for(int ix=0;ix<imgSize.width;ix++)
                {
                    iH=ix/interH;
                    if(iH>=nGridH)
                        iH=nGridH-1;
                    sum+=*ptr;
                    ptrfea1[iH]+=*ptr++;
                }
            }
            //printf("%f",sum);getchar();
            ptrfea1=ptrfea;
            for(iV=0;iV<nGridV;iV++)
            {
                for(iH=0;iH<nGridH;iH++)
                {
                    //cout<<*ptrfea1;getchar(); 
                    *ptrfea1++/=(interV*interH);
                }
            }
            /*
            if(sum<=0)
            {
                float value=sum/nGridV/nGridH;
                for(iV=0;iV<nGridV;iV++)
                    for(iH=0;iH<nGridH;iH++)
                        *ptrfea1++=value;
            }
            else
            {
                for(iV=0;iV<nGridV;iV++)
                {
                    for(iH=0;iH<nGridH;iH++)
                    {
                        *ptrfea1++/=sum;
                    }
                }
            }
            */
            //ptrfea1=ptrfea;for(iV=0;iV<nGridV;iV++)for(iH=0;iH<nGridH;iH++)printf("%f ",*ptrfea1++);getchar();
        }//end orientation
    }//end frequency
    
    // write gabor fea to files
    FILE * fpGaborFea=fopen(gaborFilename,"wb");
    if(NULL!=fpGaborFea)
    {
        //for(int i=0;i<nGaborF*nGaborO*nGridV*nGridH;i++)printf("%f ",fea[i]);getchar();
        fwrite(fea,sizeof(float),nGaborF*nGaborO*nGridV*nGridH,fpGaborFea);
        fclose(fpGaborFea);
    }
    
    delete[]fea;
    
    ::cvReleaseImage(&ptImage);
    ::cvReleaseImage(&ptTinyImage);
    ::cvReleaseImage(&ptGray);
    ::cvReleaseImage(&ptConv);
    
    
    return 0;
}
// suggested by ynli: overlap grids (slice)
int ExtractGaborSlice(const char* imageFilename, const char* gaborFilename, Gabor* gabors
, const int nGaborF, const int nGaborO, const int nSliceR)
{
    if(imageFilename==NULL || gaborFilename==NULL || NULL==gabors || nGaborF<1 || nGaborO<1 || nSliceR<1)
        return -1;
        
    IplImage*   ptImage=::cvLoadImage(imageFilename);
    if(NULL==ptImage)
        return -2;
        
    int fea_dim1;
    int fea_dim2;
    float *fea;
    ExtractGaborSlice(ptImage, &fea, fea_dim1, fea_dim2,gabors, nGaborF, nGaborO, nSliceR);
    // save fea
    //writeMatApp4BF( gaborFilename, fea, fea_dim1, fea_dim2 );
    
    delete [] fea;
        
    ::cvReleaseImage(&ptImage);
    
    return 0;
}
int ExtractGaborSlice(IplImage* ptImage, float**fea_p, int&fea_dim1,int &fea_dim2, Gabor* gabors
, const int nGaborF, const int nGaborO, const int nSliceR, const int max_size)
{
   float *fea=NULL;
    // resize to tiny image (64,64)
    int nMaxSize = ( ptImage->height > ptImage->width ) ? ptImage->height : ptImage->width;
    CvSize imgSize;
	//if (nMaxSize <= MAX_SIZE_YNLI) {
		// no need to be scaled
	//	imgSize=::cvGetSize(ptImage);
	//} else {
		// need to be scaled
		float ratio = (float)max_size / (float)nMaxSize;
		imgSize.width  = ptImage->width  * ratio;
		imgSize.height = ptImage->height * ratio;
		
	if (imgSize.width < nSliceR * 2 + 1)
		imgSize.width = nSliceR * 2 + 1;
	if (imgSize.height < nSliceR * 2 + 1)
		imgSize.height = nSliceR * 2 + 1;

	//}
    IplImage* ptTinyImage=cvCreateImage(imgSize,8, ptImage->nChannels);
    cvResize( ptImage, ptTinyImage);
        
    // convert to gray image
    //CvSize      imgSize=::cvGetSize(ptTinyImage);
    IplImage*   ptGray=::cvCreateImage(imgSize,8,1);
    if(ptTinyImage->nChannels==3)
        cvCvtColor(ptTinyImage,ptGray,CV_RGB2GRAY);
    else
        cvCopy(ptTinyImage,ptGray);
        
    IplImage*   ptConv= cvCreateImage(imgSize,IPL_DEPTH_32F, 1);
    float *ptr_conv=NULL; 

    if(imgSize.height>=(nSliceR+nSliceR+1) && imgSize.width>=(nSliceR+nSliceR+1)){   
        // if too small
         
        fea_dim1=((imgSize.height-1)/nSliceR-1)*((imgSize.width-1)/nSliceR-1); // number of slices
        //fea_dim2=4+nGaborF*nGaborO; // for x,y,r,0
		fea_dim2=6+nGaborF*nGaborO; // for x,y,r,0,w,d
        fea=new float[fea_dim1*fea_dim2];
        float *ptr_fea=NULL;
        memset(fea,0,sizeof(float)*fea_dim1*fea_dim2);
        
        int iter=0;
        int cx,cy,ix,iy;
        
        // x,y,r,0
        ptr_fea=fea;
		float *ynli_p;
		ynli_p = fea;
		int tmp_off = nGaborF*nGaborO + 6;
		int ynli_count = 0;
        for(cy=nSliceR;cy<imgSize.height - nSliceR;cy+=nSliceR)
        {
            for(cx=nSliceR;cx<imgSize.width - nSliceR;cx+=nSliceR)
            {
               // ptr_fea[0]=cx;
               // ptr_fea[1]=cy;
               // ptr_fea[2]=nSliceR;

				ynli_p[0]=cx;
                ynli_p[1]=cy;
                ynli_p[2]=nSliceR;
		
                ynli_p[4]=imgSize.width;
                ynli_p[5]=imgSize.height;

				ynli_p += tmp_off;
				ynli_count++;
            }
        }
        
        // for each frequency
        int nG=0;// which gabor filter
        int nPixels=(nSliceR+nSliceR+1)*(nSliceR+nSliceR+1);
        float sum=0;
        float f_min=99999;
        float f_max=0;
        for(int iFre=0;iFre<nGaborF;iFre++)
        {
            // for each orientation
            for(int iOri=0;iOri<nGaborO;iOri++)
            {
                nG=iFre*nGaborO+iOri;
                gabors[nG].conv_img(ptGray, ptConv, MAG);
                /*
                cvNamedWindow("conv");
                cvShowImage("conv",ptConv);
                cvWaitKey();
                cvDestroyWindow("conv");
                */
                // get slice; calculate mean; store to fea
                iter=0;
                f_min=99999;
                f_max=0;
                for(cy=nSliceR;cy<imgSize.height-nSliceR;cy+=nSliceR)
                {
                    for(cx=nSliceR;cx<imgSize.width-nSliceR;cx+=nSliceR)
                    {
                        // calculate the mean
                        sum=0;
                        for(iy=cy-nSliceR;iy<=cy+nSliceR;iy++)
                        {
                            ptr_conv=(float*)cvPtr2D(ptConv,iy,cx-nSliceR,0);
                                //cout<<*ptr_conv;getchar();
                            for(ix=cx-nSliceR;ix<=cx+nSliceR;ix++)
                            {
                                //if(*ptr_conv<1)
                                    sum+=*ptr_conv;
                                //else
                                    //sum++;
                                ptr_conv++;
                            }
                        }
                        //int temp=iter*fea_dim2+nG+4;
						int temp=iter*fea_dim2+nG+6;
                        fea[temp]=sum/nPixels;///max_mag[nG]
                                if(fea[temp]>f_max)f_max=fea[temp];
                                if(fea[temp]<f_min)f_min=fea[temp];
                        //cout<<sum/nPixels;getchar(); 
                        
                        // serrie
                        iter++;
                    }
                }/**/
                for(iter=0;iter<fea_dim1;iter++)
                {
                    //if(f_max==f_min)fea[iter*fea_dim2+nG+4]=1;
					if(f_max==f_min)fea[iter*fea_dim2+nG+6]=1;
                    //else fea[iter*fea_dim2+nG+4]=(fea[iter*fea_dim2+nG+4]-f_min)/(f_max-f_min);
					else fea[iter*fea_dim2+nG+6]=(fea[iter*fea_dim2+nG+6]-f_min)/(f_max-f_min);
                }
            }// orientation
        }// frequency
        /*
        FILE*fp=fopen("maxK.txt","w");
        for(int j=0;j<fea_dim1;j++){
            for(int i=0;i<nGaborF*nGaborO;i++)
                fprintf(fp,"%f,",fea[j*fea_dim2+i+4]);
            fprintf(fp,"\n");
        }
        fclose(fp); 
        */
    }// too small an image
    
    ::cvReleaseImage(&ptTinyImage);
    ::cvReleaseImage(&ptGray);
    ::cvReleaseImage(&ptConv);
    
    *fea_p=fea;
    fea=NULL;
    
    return 0;
}
// interface for ynli: overlap grids (slice)
int ExtractGaborSlice(IplImage* ptImage, float***fea, int&dim1,int &dim2, Gabor* gabors, const int nGaborF, const int nGaborO, const int nSliceR)
{
    if(ptImage==NULL ||  NULL==gabors || nGaborF<1 || nGaborO<1 || nSliceR<1)
        return -1;
        
    float*fea1=NULL;
	float*fea2=NULL;
	float*fea3=NULL;

	int d1, d2, d3;
	int dd1, dd2, dd3;
    
	
	ExtractGaborSlice(ptImage, &fea1, d1, dd1,gabors, nGaborF, nGaborO, nSliceR, 120);
	ExtractGaborSlice(ptImage, &fea2, d2, dd2, gabors, nGaborF, nGaborO, nSliceR, 160);
	ExtractGaborSlice(ptImage, &fea3, d3, dd3, gabors, nGaborF, nGaborO, nSliceR, 200);
	
	int nMaxSize = ( ptImage->height > ptImage->width ) ? ptImage->height : ptImage->width;

	dim1 = d1 + d2 + d3;
	dim2 = dd1;
    *fea=new float* [dim1];
    float*ptr_v_fea=fea1;
	int i, j;
	float ratio = (float)nMaxSize / 120.0;
    for(i=0;i<d1;i++)
    {
        (*fea)[i]=new float [dim2];
        for(j=0;j<dim2;j++)
        {    
			(*fea)[i][j]=*ptr_v_fea++;
		}
		for (j = 0; j < 6; j++) {
			(*fea)[i][j] = (*fea)[i][j] * ratio;
		}
    }
	delete[]fea1;

	ratio = (float)nMaxSize / 160.0;
	ptr_v_fea = fea2;
	for (i = d1; i < d1 + d2; i++) {
		(*fea)[i]=new float [dim2];
		for(j=0;j < dim2;j++) {    
			(*fea)[i][j]=*ptr_v_fea++;
		}
		for (j = 0; j < 6; j++) {
			(*fea)[i][j] = (*fea)[i][j] * ratio;
		}
	}
	delete[]fea2;
	
	ratio = (float)nMaxSize / 200.0;
	ptr_v_fea = fea3;
	for (i = d1 + d2; i < dim1; i++) {
		(*fea)[i]=new float [dim2];
		for(j=0;j < dim2;j++) {    
			(*fea)[i][j]=*ptr_v_fea++;
		}
		for (j = 0; j < 6; j++) {
			(*fea)[i][j] = (*fea)[i][j] * ratio;
		}
	}
	delete[]fea3;

    return 0;
}

// interface for ynli: overlap grids (slice)
int ExtractColorGaborSlice(IplImage* ptImage, float***fea, int&dim1,int &dim2, Gabor* gabors, const int nGaborF, const int nGaborO, const int nSliceR)
{
    if(ptImage==NULL ||  NULL==gabors || nGaborF<1 || nGaborO<1 || nSliceR<1)
        return -1;
    
	IplImage * src = cvCreateImage(cvGetSize(ptImage), 8, 3);
	cvCopy(ptImage, src, 0);
	cvCvtColor(src, src, CV_BGR2Lab);
	int nMaxSize = ( ptImage->height > ptImage->width ) ? ptImage->height : ptImage->width;

	IplImage * src3[3];
	uchar * uch[3];
	uchar *puch;
	int i, j, k, l;
	for (i = 0; i < 3; i++) {
		src3[i] = cvCreateImage(cvGetSize(ptImage), 8, 1);
	}
	for (i = 0; i < src->height; i++) {
		puch = cvPtr2D(src, i, 0);
		for (k = 0; k< 3; k++) {
			uch[k] = cvPtr2D(src3[k], i, 0);
		}
		for (j = 0; j < src->width; j++) {
			for (k = 0; k < 3; k++) {
				*(uch[k]) = *puch;				
				puch++;
				(uch[k])++;
			}
		}
	}
	
	int sizes[] = {120, 160, 200};
	//int sizes[] = {80, 120, 160};
	//int sizes[] = {120, 120, 120};
	//int nRs [] = {6, 8, 12};

	float *fea3[] = {NULL, NULL, NULL};

	vector<float *> rstFea;
	float *f_p, *f_p1;
	rstFea.reserve(1000);
	dim2 = 32 * 3 + 6;
	dim1 = 0;
	int d1, d2, d3;
	int dd1, dd2, dd3, off;
	float ratio;
	for (l = 0; l < 3; l++) {
		ratio = (float)nMaxSize / sizes[l];
		ExtractGaborSlice(src3[0], &(fea3[0]), d1, dd1, gabors, nGaborF, nGaborO, nSliceR, sizes[l]);
		ExtractGaborSlice(src3[1], &(fea3[1]), d2, dd2, gabors, nGaborF, nGaborO, nSliceR, sizes[l]);
		ExtractGaborSlice(src3[2], &(fea3[2]), d3, dd3, gabors, nGaborF, nGaborO, nSliceR, sizes[l]);
		dim1 += d1;
		
		for (i = 0; i < d1; i++) {
			f_p = new float[36 * 3 + 6];
			off = 0;
			for (k = 0; k < 3; k++) {
				f_p1 = &(fea3[k][i * dd1]);
				for (j = 0; j < dd1; j++) {
					//printf("%f ", *f_p1);
					if (j < 6 && k != 0) {
						f_p1++;
						continue;
					}
					if (j < 6) {
						f_p[off] = *f_p1 * ratio;
					} else {
						f_p[off] = *f_p1;
					}
					f_p1++;
					off++;
				}
			}
			rstFea.push_back(f_p);
			/*
			for (j = 0; j < dim2; j++) {
				printf("%f ", f_p[j]);
			}
			int  dko = 0;
			dko++;
			*/
		}
		
		for (i = 0; i < 3; i++) {
			delete [](fea3[i]);
		}
	}
	for (i = 0; i < 3; i++) {
		cvReleaseImage(&(src3[i]));
	}
	cvReleaseImage(&src);

	*fea=new float* [dim1];
	for (i = 0; i < dim1; i++) {
		(*fea)[i] = rstFea[i];
	}
	rstFea.clear();
    return 0;
}
